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Abstract — The famous Johnson-Nyquist formula relating noise 
current to conductance has a microscopic generalization re- 
lating noise current density to microscopic conductivity, with 
corollary relations governing noise in the components of the 
electromagnetic fields. These relations, known collectively in 
physics as fluctuation-dissipation relations, form the basis of 
the modern understanding of fluctuation-induced phenomena, 
a field of burgeoning importance in experimental physics and 
nanotechnology. In this review, we survey recent progress in 
computational techniques for modeling fluctuation-induced phe- 
nomena, focusing on two cases of particular interest: near-field 
radiative heat transfer and Casimir forces. In each case we review 
the basic physics of the phenomenon, discuss semi-analytical 
and numerical algorithms for theoretical analysis, and present 
recent predictions for novel phenomena in complex material and 
geometric configurations. 

Index Terms — Johnson, Nyquist, noise, fluctuation, radiation, 
heat transfer, Casimir effect, finite-difference, boundary-element, 
modeling, simulation, CAD 
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I. Introduction 

VERY electrical engineer knows the famous Johnson- 
Nyquist formula for the noise current through a resistor, 

{l^)=4kTGAf (1) 

where (P) is the mean-square noise current (Fig. la), kT 
is the temperature in energy units, G = is the conduc- 
tance of the resistor, and A/ is the measurement bandwidth. 
Equation ([T]) — which allows designers to quantify, and thus 
compensate for, the unavoidable presence of noise in physical 
circuits — is a crucial tool in the circuit designer's kit and 
a mainstay of the electrical engineering curriculum from its 
earliest stages 

Perhaps less well-known in the EE community is that 
equation ([T]) is only one manifestation of a profound and 
far-reaching principle of physics — the fluctuation-dissipation 
theorem — that relates the mean-square values of various fluc- 
tuating quantities to certain physical parameters (known as 
generalized susceptibilities) associated with the underlying 
system. In equation ([T]), the fluctuating quantity is the noise 
current through the resistor, and the generalized susceptibility 
is the conductance; more generally, as we will see below, the 
fluctuation-dissipation concept allows us to quantify fluctu- 
ations not only in macroscopic device currents but also in 
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microscopic current densities, from which it is a short step 
to obtain fluctuations in the components of the electric and 
magnetic fields inside and outside material bodies (Fig. lb). In 
this case, we will see that the key tools turn out to be nothing 
but the familiar dyadic Green 's functions, which describe the 
electromagnetic fields of prescribed current sources and are 
computable by any number of standard methods of classical 
electromagnetism. It is remarkable that many problems in the 
field of fluctuation-induced phenomena, which would at first 
blush seem to necessitate complex statistical-mechanical and 
quantum-mechanical reasoning, in fact reduce in practice to 
applications of classical electromagnetic theory that would be 
familiar to any electrical engineer. 

But why would we want to quantify noise in the individual 
components of the electromagnetic fields around material bod- 
ies? The answer is that these microscopic field fluctuations can 
mediate macroscopic transfers of energy or momentum among 
the bodies, which become especially dramatic for bodies at 
submicron separations. In the former phenomenon — near-field 
radiative heat transfer — fluctuating fields in micron- scale gaps 
between inequal-temperature bodies can lead to a rate of heat 
transfer between the bodies that can drastically exceed the rate 
observed at larger separations |2|. In the latter phenomenon — 
the Casimir effect — fluctuating fields around bodies give rise 
to attractive and repulsive forces between the bodies, which 
generalize the familiar van der Waals interactions between 
molecules ||3|. Both phenomena become negligibly small for 
bodies separated by distances of more than a few microns, 
which places their observation squarely within the domain of 



(a) 




Fig. 1. From macroscopic to microscopic noise, (a) The current through a 
resistor exhibits thermal noise with mean-square amplitude proportional to the 
conductance [the Johnson-Nyquist formula, equation ([^]. (b) More generally, 
the microscopic current density inside a slab of conducting material exhibits 
fluctuations with mean-square amplitude proportional to the microscopic 
conductivity [the fluctuation-dissipation theorem, equation |2j]. Knowledge 
of these microscopic current fluctuations, together with the dyadic Green's 
functions of the system, allow us to predict the mean-square fluctuations in 
the components of the electromagnetic fields in space [equations ^ and fT2)]. 
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Fig. 2. A selective timeline indicating the most complex geometries for which rigorous calculations of Casimir interactions (upper) or near-field radiative 
heat transfer (lower) were possible at various historical epochs. Note that computational techniques such as finite-difference grids and boundary-element 
discretization, which have been used in electrical engineering for decades, have only been introduced to the study of fluctuation-induced phenomena within 
the past five years. 



nanoscale physics and engineering. 

Although the study of electromagnetic field fluctuations has 
been an active area of physics for decades, its relevance to 
electrical engineering was limited for most of that time to 
equation ([T]) and other relations quantifying noise in circuits. 
In the past fifteen years, however, this situation has begun to 
change; advances in fabrication and measurement technology 
have ushered in a golden age of experimental studies of 
fluctuation-induced phenomena |j2|, p7| , and there is reason 
to believe that this fledgling field of experimental physics 
will soon become relevant to electrical engineering in areas 
such as thermal lithography and MEMS technology. This 
experimental progress has created a demand for modeling and 
simulation tools capable of predicting fluctuation phenomena 
in realistic experimental configurations, including the complex, 
asymmetric geometries and imperfect materials present in real- 
world systems. 

The evolution of theoretical tools for modeling fluctuation- 
induced phenomena mirrors the historical development of 
techniques for solving classical electromagnetic scattering 
problems. In the latter case, the earliest calculations were 
restricted to highly symmetric geometries (such as Mie's 1908 
treatment of scattering from spheres) for which a conve- 
nient choice of coordinates and special-function solutions of 
Maxwell's equations allow the problem to be solved ana- 
lytically (or at least semi- analytically — that is, with results 
obtained as expansions in special functions, which in practice 
are then evaluated numerically |18|). Later, fully numerical 
techniques capable of handling more general geometries grad- 
ually became available, including the finite-difference, finite- 
element, and boundary-element methods introduced in the 
1960s, and today the problem of electromagnetic scattering is 
addressed by a wealth of comprehensive off-the-shelf CAD 
tools capable of handling extremely complex material and 
geometric configurations. 

Advances in the modeling of near-field radiative transfer 
and Casimir phenomena have proceeded in similar order (Fig. 
|2]). In both cases the first calculations were restricted to the 
simplest parallel-plate geometries |^, | |T2| , |[T3|; these were 



later extended to other simple shapes such as cylinders 1 8 1 and 
spheres f5l, |9l, |T4| , |T9|-pTj, and, more recently, tools for 
general geometries have become available p2|-p4|. All of 
these developments, however, have lagged their antecedents 
in the classical-scattering domain by many decades; indeed, 
even for the relatively simple case of two interacting spheres, 
the Casimir force was only calculated in 2007 |9| and the 
near-field radiative transfer only in 2008 1 14|. One reason for 
this lag is the relative paucity of experimental data, which — 
as noted above — are significantly more difficult to gather for 
fluctuation-induced phenomena than for classical scattering. 
But perhaps the main reason that practical computations of 
fluctuation-induced phenomena have been so long in coming is 
simply that the problems present extraordinary computational 
challenges. Indeed, as we will see below, calculations of 
near-field radiative-transfer and Casimir phenomena may be 
reduced in practice to the solution of classical scattering 
problems — but a great number, thousands or even millions, 
of separate scattering problems must be solved to compute 
the heat transfer or Casimir force for a single geometric 
configuration. 

As a result, algorithms for predicting fluctuation phenomena 
tend to start with techniques familiar to electrical engineers — 
including the T-matrix, finite-difference, and boundary- 
element methods of computational electromagnetism — but 
then proceed to combine and modify these techniques in novel 
ways to obtain computational procedures that can run in a 
reasonable length of time. The goal of this review is to describe 
these computational techniques — and some of the results that 
they have predicted — in ways that will make sense to electrical 
engineers. 
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II. Fluctuations in Electromagnetic Sources and 
Fields: The Johnson-Nyquist Law and Beyond 

The microscopic generalization of equation ([T]) is fTSj , p5| 

yi(a;,x)J*(cj,xO^ 



1 



— coth —-- 
2 \2kT) 



Cr(cJ,x) (2) 



where Ji(co',x) is the i\h cartesian component of the 
microscopic current density at position x and frequency 
cj, ?i is Planck's constant, and cr(a;,x) is the position- 
and frequency-dependent conductivity, [a is related to 
the imaginary part of the dielectric permittivity according 
to cr (cj, x) = a; • Im e(a;, x); here and throughout we assume 
that e is linear and isotropic] In signal-processing language 
familiar to electrical engineers, the right-hand side of equation 
^ is the power spectral density (PSD) of a colored-noise pro- 
cess; the fact that the PSD is frequency-dependent (i.e. the fact 
that this is "colored" instead of white noise) corresponds, in 
the time domain, to the nonvanishing of correlations between 
currents at nearby time points. 

The similarity between equations ([T]) and ^ is obvious: on 
the left-hand side we have a mean product of currents, while 
on the right-hand side we have a temperature-dependent factor 
and a measure of conductivity. However, the microscopic 
equation ^ extends the macroscopic equation ([T]) in two 
important ways. 

First, whereas equation ([T]) is a low-frequency, high- 
temperature approximation that neglects quantum-mechanical 
effects, equation ^ remains valid at all temperatures and 
frequencies and explicitly includes quantum-mechanical ef- 
fects. Indeed, taking the low-temperature limit of the bracketed 
factor in ([2]), we find 

Tiuj 



I'll 11 



Tiuj f Tiuo \ 
coth -—- 
\2kT J 



(3a) 



and equation ^ thus predicts non-zero current fluctuations 
even at zero temperature: the well-known quantum-mechanical 
zero-point fluctuations. In the high-temperature limit, on the 
other hand, we have 



lim 

T^oo 



huj f Tiuj \ 
— coth -—- 
2 \2kT) 



(3b) 



this is the classical regime, in which all dependence on ft is 
lost and we recover the simple linear temperature dependence 
of ([T]). The classical regime is defined by the condition 

a requirement that in practice is always satisfied in circuit- 
design problems, but which may be readily violated for 
infrared and optical frequencies {uj > 10^^ rad/sec). 

The second way in which equation ^ extends the reach of 
the Johnson-Nyquist result is that, whereas ([T]) describes only 
macroscopic currents, ^ gives information on the microscopic 
current density, which in turn can be used to predict fluctua- 
tions in the components of the electric and magnetic fields. The 
relevant tools for this purpose are the dyadic Green's functions 
(DGFs), well-known to electrical engineers from problems 



ranging from radar and antenna design to microwave device 
modeling |18|. To recall the definition of these quantities, 
suppose we have a material configuration characterized by 
spatially-varying linear permittivity and permeability functions 
{e(co', x), /i(cj, x)}. (In most of the problems we consider, e 
and /i will be piecewise constant in space.) Then the electric 
DGF describes the field due to a point source in the presence 
of the material configuration: 

Gf^(e,/i;a;;x,x') 

_ /z-component of electric field at x due to\ 

~ \3. j -directed point current source at x' y ^ ^ 

while the magnetic DGF similarly gives the magnetic 
field of a point current source. (Here and throughout, all 
fields and currents are understood to have time dependence 
~ e~*^^) In ^ we have indicated the dependence of G on 
the spatially- varying material properties e and /i; the DGFs for 
a given material configuration can be computed using standard 
techniques in computational electromagnetism, after which the 
fields at arbitrary points in space due to a prescribed current 
distribution may be computed according to 

Ei{uj,x.) = j G%{uj]yi,y:!)Jj{uj,^)d-x! (6a) 

Hi{uj,^) = j Gfj{uj',^,^)Jj{uj,^)d^. (6b) 

Note that the long-range nature of the G dyadics ensures that 
the fields are nonvanishing even at points x in empty space, 
i.e. points at which there are no currents or materials. 

Armed with equations ^ and ([6]), we can now make 
predictions for noise in the components of the electromagnetic 
fields. For example, the mean Poynting flux at a point x is a 
sum of terms of the form (with uj arguments to E, G, and J 
suppressed) 

(Ei(x)F;(x)) 

= j dx'dx"Gffe(x,x')G,t(x,x")(Jfc(x')JKx")) 
Inserting ([2]), 

d^GU^, x')Gf,*(x, x')e [uj, T(x')]ct(w, x') (7) 



/ 



where T(x) is the local temperature and 6[cj,T] = 
^ coihTiuj / 2kT is the statistical factor in equation 
(Summation over repeated tensor indices is implied here and 
throughout.) 

The obvious advantage of an equation like ^ is that it 
reduces a problem in quantum statistical mechanics (deter- 
mination of the electromagnetic field fluctuations at x) to a 
problem in classical electromagnetic scattering (computation 
of the DGFs G^'^). The difficulty of this approach lies in 
the great number of scattering problems that must be solved. 
Indeed, equation ^ says that, to compute the Poynting flux 
at a single point x, we need the DGFs connecting x to all 
points x' at which the conductivity is nonvanishing; for a 
typical problem involving two dissipative bodies in vacuum, 
this amounts to a solving a separate scattering problem for 
each point in the volume of each body. Moreover, even 
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after completing all of these calculations we have still only 
computed the Poynting flux at a single point x; in general 
we will want to integrate this flux over a surface to get the 
total power transfer at a given frequency, and subsequently to 
integrate over all frequencies to get the total power transfer. 

Thus the fluctuation-dissipation concept, in the form of 
equations ([2]) or (|7]), performs the great conceptual service 
of reducing predictions of noise phenomena to problems of 
classical electromagnetic scattering, but leaves in its wake the 
practical problem of how to solve the formidable number 
of scattering problems that result. This difficulty has been 
addressed in a variety of ways, some of which we will review 
in the following sections. 

III. Near-Field Heat Radiation: 
Fluctuation-Induced Energy Transfer in 
Nanoscopic Systems 

Fluctuating currents in finite-temperature bodies give rise 
to radiated fields which carry away energy. If there are other 
bodies (or an embedding environment) present at the same 
temperature, then any energy lost by one body to radiation 
is replenished by an equal energy absorbed from the radia- 
tion of other bodies. However, between objects at different 
temperatures there is a net transfer of power, whose rate we 
can calculate in terms of the temperatures and electromagnetic 
properties of the bodies. 

Historically, the first step in this direction was the Stefan- 
Boltzmann law, a triumph of 19th-century physics which held 
that the power radiated per unit surface area of a temperature- 
T body was simply t/ctsbT^, where (Jsb is a universal constant 
and T], the emissivity, is a dimensionless number between 
and 1 characterizing the electrical properties of the body 
(specifically, its propensity to emit radiation relative to that 
of a perfect emitter or black body). The Stefan-Boltzmann 
prediction is based on an approximation that simplifies the 
electromagnetic analysis: it considers only propagating elec- 
tromagnetic waves, neglecting the evanescent portions of the 
E and H fields that exist in the vicinity of object surfaces. This 
is a good approximation when computing the power transfer 
between a single body and its environment, or between two 
inequal-temperature bodies separated by large distances. 

However, when inequal-temperature bodies are separated by 
short distances, evanescent fields can contribute significantly 
to the Poynting flux and the rate of power transfer may deviate 
significantly from the Stefan-Boltzmann prediction. The length 
scale below which distances are to be considered "short" is the 
thermal wavelength, 

. ?ic /300 K\ 

Ar=^«7.6Mm-^— j, 

and thus, in practice, observing deviations from the Stefan- 
Boltzmann law requires measuring the heat flux between two 
bodies maintained at inequal temperatures and at a surface- 
surface separation of a few microns. This formidable experi- 
mental challenge has recently been met by several groups | 2J, 
| [26| , and this progress has spurred the development of new 
theoretical techniques for predicting the heat flux between 



closely-spaced bodies with realistic material properties and 
various shapes, which we now describe. 

A. Radiative Heat Transfer as a Scattering Problem 

Consider two homogeneous bodies Si, 2 separated by a 
short distance and maintained at separate internal thermal 
equilibria at temperatures Ti 2. (We will consider the bodies 
to exist in vacuum; the case of a finite-temperature embedding 
environment is a straightforward generalization.) The rate at 
which energy is absorbed or lost by body 1 is given as a 
surface integral of the mean Poynting flux, 

Pi{uj) = \j^ (e(c^,x) xH*(a;,x))-dS, (8) 

where Si is the surface of body 1 (or, equivalently, a fictitious 
bounding surface containing body 1 and no other bodies) and 
dS is the inward-pointing surface normal. Applying equa- 
tion ([7]) reduces the quantity in brackets to integrals over the 
volumes of the bodies (again suppressing uj arguments to G): 

Pi{^) (9) 

+ a2(c^)e[c^,T2] ^ Gf^(x,xOGf/(x,x0^ix l^/^fc. 

where cri^2 are the conductivities of the bodies. [Here we have 
used the Levi-Civita symbol £ijk to write the components 
of the cross product as (A x B)/e = SijkAiBj.] Note that 
equation ^ includes integrations over the volumes of both 
bodies, since there are fluctuating sources present in both 
bodies. Intuitively one might expect that reciprocity arguments 
could be exploited to relate the two terms to one another and 
hence streamline the calculation to involve integration over 
just one body; this intuition is indeed born out in practice, as 
discussed below |15|. 

Equation ^ reduces the calculation of the net energy 
transfer to or from a body to the classical electromagnetic 
scattering problem of computing the DGFs for a geometry 
consisting of our two material bodies Si, 2- The difficulty, as 
anticipated above, is that we must solve a great number of 
scattering problems; in principle, for each surface point x 
and each volume point x' in the combined surface-volume 
integrals in ^ we must solve a separate scattering problem 
(computing the fields at x due to individual point sources at 
xO- This challenge is in fact so formidable that computations 
for geometries even as simple as two spheres have only 
become available in the past few years, using techniques which 
we now review. 

B. Semi- Analytical Approaches to Radiative Transfer 

A first strategy for evaluating ^ is to consider certain 
highly symmetric geometries for which a convenient choice of 
coordinates allows the DGFs to be evaluated analytically. For 
example, the earliest near-field heat-transfer calculations 1 12], 
|[T3l took the two objects to be semi-infinite planar slabs, 
in which case the DGFs are analytically calculable. More 
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recently, several groups have extended this approach to other 
highly symmetric geometries in which special-function ex- 
pansions of the DGFs are available |14|, |27|-|32|. A par- 
ticularly convenient tool here is the "matrix" approach to 
electromagnetic scattering, a technique first discussed in these 
PROCEEDINGS in 1965 |33|. To solve scattering problems 
in this approach, one begins by writing down two sets of 
functions, {E5^(x)} and {E^^^(x)}; these are solutions of 
Maxwell's equations, in an appropriate coordinate system, 
which respectively describe electromagnetic waves propagat- 
ing inward from infinity to our scattering geometry and 
outward from the scatterer into open space. (For example, 
in spherical geometries the {E^} will be products of vector 
spherical harmonics and spherical Bessel functions fTS].) The 
disturbance in the electromagnetic field due to a scattering 
object is then entirely encapsulated in the object's T-matrix, 
denoted T, whose m, n element gives the amplitude of the mth 
outgoing wave for a scatterer illuminated by the nth incoming 
wave. In other words, 

if the incident field is E^^^(x) = Ej^(x) 
then the scattered field is E'"^^(x) = ^T^nE^'(x). 




frequency co (liic/ a) 

Fig. 3. Near- field radiative heat transfer between patterned and unpatterned 
SiC slabs |15|. The solid black curve plots the spectral density of power 
flux between SiC photonic crystals (inset) maintained at inequal temperatures 
and surface-surface separation d. The dashed red curve plots the power flux 
between unpatterned SiC slabs. (In both cases, the power flux is normalized 
by the power flux that would obtain between the same structures at infinite 
separations (i — ^ oo.) 



Because the T-matrix for a body encodes all information 
needed to understand its scattering properties, it is often 
possible to express the solution to a radiative-transfer problem 
in terms of simple matrix operations on the T-matrices of 
the objects involved. As an illustration of the sort of concise 
expression that can result from this procedure, the methods 
of Ref. |27| lead to a simple trace formula for the spectral 
density of heat radiation from a single sphere at temperature 
T (341: 



H{uj, T) = -26' [cj, T] ^ I Re 

n ^ 



(10) 

where T is the T-matrix of the sphere, the sum runs over its 
diagonal elements, and 6' is just 9 minus the contribution of 
the zero-point energy term. 



The obvious advantage of an equation like ([TO]) is that 
it is simple enough to be implemented in a few lines of 
MATHEMATICA or MATLAB for objects whose T-matrix is 
known analytically. The difficulty is that there are not many 
such objects; indeed, the only lossy scatterers for which the 
T-matrix may be obtained in closed form are spheres, infinite- 
length cylinders, and semi-infinite slabs. (Idealizing the mate- 
rials as lossless metals extends the list of shapes for which the 
T-matrix is known analytically |35J, but this is not useful for 
radiative-transfer problems because lossless materials neither 
absorb nor radiate energy.) To make predictions for shapes 
outside this narrow catalog we must turn instead to numerical 
methods. 

C. Numerical Approaches to Radiative Transfer 

One approach to numerical heat-transfer modeling is to 



combine matrix-trace formulas in the spirit of equation (10) 
with a numerical technique for computing the T-matrices 



of irregularly- shaped objects. This technique was pursued in 
Ref. fT6| , which investigated heat transfer from hot tips of 
various shapes to a cool planar substrate at micron-scale dis- 
tances. In this work, a boundary-element scattering code was 
used to compute numerical T-matrices for finite cylinders and 
finite-length cones; a surprising conclusion was that conical 
tips, despite tapering to a point, nonetheless exhibit less spatial 
concentration (i.e. a larger and more diffuse spot size) of heat 
transferred to the substrate as compared to cylindrical tips. 

An alternative numerical approach to heat-transfer calcu- 
lations is to bypass the T-matrix approach in favor of a 
more direct assault on equation ([9]) fT5| , | [36| ; here a "brute- 
force" approach can deliver great generality with minimal 
programming time, at the expense of much computer time. 
Physically, the situation described by equation ([9]) is that 
we have randomly fluctuating currents distributed throughout 
the interior of our material bodies, and we wish to compute 
the fields to which these currents give rise. A particularly 
convenient way to do this computation is to run a time-domain 
simulation, in which we calculate the fields due to a random 
time- varying current distribution whose correlation function in 
the frequency domain satifies equation ([2]); by repeating this 
calculation for many randomly-generated current distributions 
and averaging the results, we obtain approximate ensemble 
averages of the time-domain E and H fields, which we 
may then Fourier- analyze to obtain frequency spectra. This 
approach is rendered computationally feasible by exploiting 
several properties of equation ^ and of Maxwell's equations. 
First, absence of spatial correlation: the 6 function in ^ en- 
sures that currents at different locations in space (in particular, 
currents in different bodies) are uncorrelated and may thus be 
chosen to have independent random phases. Second, linearity: 
although equation ^ calls for stochastic currents with non-flat 
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spectral density shaped by the factor 6[a;, T] — what engineers 
might think of as "colored noise" — the linearity of Maxwell's 
equations ensures that we can instead compute the fields due to 
white-noise currents, which are significantly easier to generate 
in the time domain, and only later multiply the resulting 
frequency spectrum by the appropriate shaping factor. Finally, 
reciprocity: the flux absorbed by body B2 due to radiating 
sources in Bi is equal to the flux absorbed by Bi due to sources 
in S2. This observation allows us to place our stochastic 
sources only in Bi and compute the resulting flux only into 

Combining these arguments leads to a simple expression for 
the spectral density of the net heat flux between bodies p5| : 

H{uj, Ti , T2) = ^{uj)[q[uj, Ti] - e [a;, T2] } 

where ^ is the flux into one of the objects due to random 
(white-noise) current sources in the other object. In practice, ^ 
is computed using a finite-difference time-domain technique, 
with random current sources placed at grid points throughout 
the volume of the bodies and the results averaged over many 
(~ 60) simulations. 

Fig. |3] illustrates the type of result that may be obtained 
using this method 1 15 |. The solid curve in the figure plots the 
spectral density of power flux between two one-dimensional 
photonic crystals of SiC separated by a short distance d (inset). 
The dashed curve plots the power flux between unpatterned 
SiC slabs. (In both cases, the power flux is normalized by the 
flux between the same structures at infinite separation d 
00.) The patterning of the slabs drastically modifies the flux 
spectrum as compared to the unpatterned case. 



IV. Casimir Forces: Fluctuation-Induced 
Momentum Transfer in Nanoscopic Systems 

In the previous section, we considered applications of 
fluctuation-dissipation ideas to situations out of thermal equi- 
librium, and we noted the fierce computational challenges 
that arise from the need to solve separate scattering problems 
for each point in the volume integration in (|7]). At thermal 
equilibrium, a major simplification occurs which significantly 
reduces computational requirements. The situation is most 
clearly displayed by considering the mean product of E-field 
components, which reads, in close analogy to (|7]), 



(i?.(x)£:;(x')) 



(11) 



J dy Gffc(x, y)G%{^', y)e [to, T{y)] <t(w, y) 



The key point is that, at thermal equilibrium, T(y) = T 
is spatially constant, whereupon the statistical factor may be 
pulled out of the integral to yield 



6[c^,T] I ^yGf,(x,y)Gf:(x^y)a(a;,y) 



= -6[c^,T]Im Gf,(x,xO. 



(12) 



equations p7|.) Thus, evaluating a mean product of field 
components at thermal equilibrium requires the solution of 
only a single scattering problem, in contrast to the formally 
infinite number of scattering problems required for out-of- 
equilibrium situations. 

Of course, the heat-transfer calculations of the previous 
section are not very interesting at thermal equilibrium, in 
which by definition there can be no net transfer of energy 
between bodies. However, a different sort of fluctuation- 
induced phenomenon — the Casimir effect — gives rise to non- 
trivial interactions among bodies even at the same temperature 
(and even at zero temperature), and constitutes a second major 
branch of the study of electromagnetic fluctuations. 



A. The Casimir Effect 

In 1948 jSSj, Casimir and Polder generalized the van der 
Waals (or "London dispersion") force between fluctuating 
dipoles of molecules and other small particles, which depends 
on the distance r between the particles like l/r^, to a "re- 
tarded" force that varies like at large distances (typically 
tens of nanometers) where the finite speed of light must be 
taken into account. Later that year [fl, Casimir considered the 
region between two parallel mirrors as a type of electromag- 
netic cavity, characterized by a set of cavity-mode frequencies 
{cjo{d)} depending on the mirror separation distance d. By 
summing the zero-point energies [equation ([3^)] of all modes 
and differentiating with respect to d, Casimir predicted an 
attractive pressure between the plates of magnitude 

F Tr'^hc 10-^ atm 

(13) 
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negligible at macroscopic distances but significant for surface- 
surface separations below a few hundred nanometers. 

The Casimir effect was subsequently reinterpreted |[39j, 
I4Q} as an interaction among fluctuating charges and currents 
in material bodies, a perspective which allows the use of 



(In going to the last line here we used a standard identity in 
electromagnetic theory which follows directly from Maxwell's 



fluctuation-dissipation formulas like ([12]) to predict Casimir 
forces in situations where the cavity-mode picture would be 
unwieldy. In fact, the Casimir effect has been interpreted in 
a bewildering variety of ways; in addition to the zero-point- 
energy picture of Ref. |4 1 and the source- fluctuation picture of 
Ref. p9| , there are path-integral formulations | |23| , world-line 
methods f4T\, and ray-optics approaches |42l, to name but a 
few. Each of these perspectives emphasizes different aspects 
of the underlying physics, although of course all physical 
interpretations lead ultimately to mathematically equivalent 
final results |24|. However, despite the plethora of theoretical 
perspectives, and even with the simplifications afforded by 
thermal equilibrium, the calculations remained so challenging 
that force predictions for all but the simplest geometries were 
practically out of reach, and — with experimental progress 
hampered by the difficulty of measuring nanonewton forces 
between bodies at sub-micron distances — for many decades 
there was little demand for computational Casimir methods 
that could handle general geometries and materials. 

This situation began to change about 15 years ago with the 
advent of precision Casimir metrology |[43|, and since that 
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time the Casimir effect has been experimentally observed in 
an increasingly wide variety of geometric and material config- 
urations (for recent reviews of experimental Casimir physics, 
see fT?! , |[44|). This experimental progress has spurred the 
development of theoretical techniques capable of predicting 
Casimir forces and torques in complex, asymmetric geometries 
with realistic materials, which we now review. 

B. The Casimir Effect as a Scattering Problem 



As in Section III- A we consider two bodies Si, 2 in vacuum. 
In equation ([8| we integrated the average Poynting flux over a 
surface surrounding a body to obtain the rate of energy transfer 
to that body. To compute the rate of momentum transfer to the 
body — that is, the force on the body — we proceed analogously, 
but now instead of the Poynting flux we integrate the average 
Maxwell stress tensor: 

T{uj) = J^(T{uj,^)ydS, (14) 

where the components of T are given in terms of the compo- 
nents of E and H as 



2 



Inserting ^V2\ and its magnetic analogue into ([14]) now yields 
an expression analogous to ([9]) — but simplified by the absence 
of volume integrals — which at temperature T = takes the 
form, for the i component of the force, 

= ^Im^ |eoaf^(cJ,x,x) +/ioSl^(c^,x,x) (15) 
Gkk (^^ + Mo Gkk (^^ 



2 



Here ^(x, x'), the scattering part of a DGF G, is the con- 
tribution to G which remains finite as x' ^ x; this is just 
the field at x due to currents induced by a point source at x^ 
but neglecting the direct contribution of that point source. [In 
( 15 ), is the scattering part of the DGF that relates magnetic 
fields to magnetic currents.] 

Equation ( 15 ), like equation ([9|, reduces our problem to that 
of determining the DGFs for our material configuration, and in 
principle we could now proceed to evaluate the surface integral 
in ( p3] ) with the integrand computed by standard scattering 
techniques. For Casimir calculations, however, the situation is 
complicated by an important subtlety, which we now discuss. 

C. Transition to the Imaginary Frequency Axis 

In contrast to the heat-transfer problems discussed in the 
previous section, for Casimir problems we will not typically 
be interested in the contributions of individual frequencies 
but will instead seek only the total Casimir force on a body. 



(16) 



obtained by integrating ^T5\ over all frequencies: 

POO 

Fi= Ti{uj)duj. 
Jo 

But naive attempts to evaluate equation ([16]) numerically are 
doomed to failure by the existence of rapid oscillations in the 
integrand, as pictured in Fig. |4^ for the particular case of 




2 3 4 5 

Real frequency co (units of 27Tc/a) 



4Ci 



is 





/ Im cj ) 






► 

^ Re oj 
^ X 
X 


f Complex LJ plane | 

(b) 



0.5 1 1.5 2 

Imaginary frequency ^ (units of 27Tc/a) 

Fig. 4. Transition to the imaginary frequency axis, (a) As a function 
of real frequency lu, the Casimir force integrand ^{uj) of equation (Ts) — 
shown here for the case of parallel metallic plates separated by a distance a 
(inset) — exhibits severe oscillations which effectively prohibit evaluation of 
the integral ([T6j by numerical quadrature. (For the particular case of parallel 
plates, the expression for the Casimir force integrand is known as the Lifshitz 
formula |40|.) These oscillations are associated with cavity resonances, which 
show up mathematically as poles in the lower half of the complex frequency 
plane (inset); the real part of the pole corresponds to the resonance frequency, 
while the imaginary part corresponds to the width (or the inverse lifetime) 
of the resonance, (b) Rotating to the imaginary frequency axis (inset) moves 
the contour of integration away from the cavity-resonance poles, resulting in 
a smooth integrand that succumbs readily to numerical quadrature. 



the Casimir force between parallel metallic plates in vacuum. 
The origin of these oscillations is not hard to identify: they 
are related to the existence of electromagnetic resonances in 
our scattering geometry, which correspond mathematically to 
poles of the integrand in the lower half of the complex uj plane. 
(The oscillatory nature of the force spectrum was emphasized 
in Ref. (45} , and the implications for numerical computations 
were discussed in Ref. ||46|.) 



But this diagnosis of the problem suggests a cure: thinking 
of ([16]) as a contour integral in the complex frequency plane, 
we simply rotate the contour of integration 90 degrees and 
integrate over the imaginary frequency axis (Fig. 4b). This 
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procedure, known in physics as a Wick rotation | [47| , yields 

F^= / J'^{Od^ (17) 

Jo 

where uj = and T now involves the DGFs evaluated at 
imaginary frequencies: 

m) = ^ l^[eoglj{^,^,^) + f^oGfj{^,^,^) (18) 



2 



eo Sfe/c , X, x) + /io Qkk . X, x) 



The Wick rotation is possible here because the DGFs are 
analytic functions in the upper half of the complex uj plane. 
This is a well-known consequence of causality: the fields 
arise after the current fluctuations that generate them |48|. 
Another consequence of causality is that, for passive materials, 
the permittivity and permeability functions on the imaginary 
frequency axis {e(z^), guaranteed to be real-valued 

and positive | [49| . 

Physically, the transition to the imaginary frequency axis 
corresponds to replacing the oscillatory time dependence 
~ g-*u;t q£ fields and currents with an exponentially 
growing time dependence ~ e~^^^; for frequency-domain com- 
putational electromagnetism, this has the effect of replacing 
the spatially oscillatory Helmholtz kernel C ) with an 
exponentially decaying kernel (^^^)- As illustrated in Fig. 
4b, the imaginary-frequency Casimir force integrand 
is a well-behaved smooth function that succumbs readily to 
numerical quadrature. 

Equations ^Vl) and (18) are valid at zero temperature. At 
finite temperatures T > 0, we must include a factor 6[^, T] ^ 
cothih^/2kT under the integral sign; in this case, it is well- 



known in physics |40 1 that the integral ( 17 ) over the imaginary 
frequency axis may be evaluated using the method of residues 
to obtain 

2..T.,_,., (19) 



n 



n=0 



where = 2n7rkT/h, the Matsubara frequencies, are just 
the poles of the coth factor on the imaginary frequency axis. 
[The primed sum in ([19]) indicates that the n = term enters 
with weight 1/2.] Computationally, the upshot of equation ([T9| 
is that finite-temperature Casimir forces are computed with 
no more conceptual difficulty than zero-temperature forces. 



with the integral in (17) simply replaced by the sum in 



(19), although the need to evaluate equation (18) in the 
limit of zero frequency = 0+) poses challenges for some 
methods of computational electromagnetism [5QJ , |^ |. The 
temperature dependence of Casimir interactions is a topic of 
recent theoretical |52j and experimental [ [53| interest. 

D. Semi- Analytical Approaches to Casimir Computations 

Like the first studies of near-field radiative transfer, the first 
generation of theoretical Casimir techniques focused on highly 
symmetric geometries for which analytical scattering solutions 
are available |^, ||28l, ||54l-|pS|. As an example of the 
type of concise expression that may be obtained via these 



methods, the zero-temperature Casimir force between two 
compact bodies with center-center separation vector R may 
be expressed in the form |9 | 



Tr 



h 

where the matrix M has the block structure 

T-^ U(R) \ 
Ut(R) T^-^ J 



(20) 



M : 



here is the T-matrix for body n and U(R) is a translation 
matrix, which relates spherical Helmholtz solutions about dif- 
ferent origins and for which closed-form analytic expressions 
are available |59|. [The partial derivative in ( [2Q| ) is taken with 
respect to a rigid displacement of one body in the ith cartesian 
direction.] 



Like equation ([TO]), the formula ( [20] ) is simple enough that 
it can be implemented in just a few lines of MATHEMATICA 
or MATLAB code for geometries in which the T— matrix is 
known analytically. Again, however, such geometries are rare, 
and for more complicated geometric configurations we must 
turn to numerical methods. 

E. Numerical Approaches to Casimir Computations 

The most direct way to apply numerical techniques to 
Casimir computations is simply to evaluate the surface in- 
tegral in ( [T8] ) by numerical cubature, with the Q tensors 
at each integrand point x evaluated by solving a numerical 
scattering problem in which we place a point source at x 
and compute the scattered fields back at the same point x. 
In principle, this scattering problem may be solved by any 
of the myriad available techniques for numerical solution 
of scattering problems (although the need for imaginary- 
frequency calculations poses something of a limitation in 
practice). To date, computational Casimir methods based on 
numerical evaluation of ( [TSj ) have been implemented using a 
variety of standard techniques in computational electromag- 
netism: the finite-difference frequency-domain method | [46| , 
1 60 1, the finite-difference time-domain method [with some 
transformations to convert the integral over frequencies in ( [TSj ) 
into an integral over the time-domain response of a current 
pulse] |[6T|-|[63|, and the boundary-element method |[64|, (651. 



Compared to the special-function approaches discussed in 
Section [IV-D[ any one of these numerical methods offers the 
significant practical advantage of handling arbitrarily complex 
geometries with little more difficulty than simple geometries. 
Among the various numerical methods, the finite-difference 
methods have the advantage of greater generality — in the 
sense that they can readily handle arbitrarily complex material 
configurations, including anisotropic and continuously- varying 
dielectrics — while the boundary-element methods have the ad- 
vantage of greater computational efficiency for the piecewise- 
homogeneous material configurations typically encountered in 
practice. 

As an illustration of the type of problem that is facilitated by 
numerical Casimir methods. Figure [5] plots the force between 
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Fig. 5. Casimir force between elongated pistons confined between parallel 
plates |46|. The lower inset depicts the geometry, while the upper inset shows 
the finite-difference grid used to model the cross-section of this z-invariant 
structure. The force between the pistons exhibits a suprising non-monotonic 
dependence on the separation distance h between the pistons and the plates. 
[The quantity plotted is the actual force divided by the proximity-force 
approximation (PFA) to the force, a convenient /i-independent normalization.] 



elongated square pistons confined between parallel plates (all 
bodies are perfect conductors), as computed using a finite- 
difference technique |46|. The lower inset in the figure depicts 
the geometry, while the upper inset shows the finite-difference 
grid used to model the cross-section of this z-invariant struc- 
ture. The force between the pistons exhibits a suprising non- 
monotonic dependence on the separation distance h between 
the pistons and the plates. 

F. Fluctuating-Surf ace-Current Approach to Casimir Compu- 
tations 

The finite-difference and boundary-element methods de- 
scribed above have the advantage of great generality, in that 
they treat bodies of arbitrarily complex shapes with no more 
difficulty than simple symmetric bodies. However, the need 



for numerical evaluation of the surface integral in (18) adds 
a layer of conceptual and computational complexity that is 
absent from the concise expression ( [2Q| . 

An alternative is the recently &t\t\o^t& fluctuating- surface- 
current approach pQ| , |[66|-|[68|. In the FSC technique, we 
begin with a boundary-element-method (BEM) approach to 



evaluating the DGFs in ([18]). Instead of proceeding numeri- 
cally, however, we exploit the structure of the BEM technique 
to obtain compact analytical expressions for the DGFs in 
fully-factorized form, involving products of factors depending 
separately on the source and evaluation points. Inserting these 
expressions into ([18]) then turns out to allow the surface 
integral to be evaluated analytically, in closed form, leaving 
behind only straightforward matrix manipulations |[67j, | [68| . 
The final FSC formula for the Casimir force. 



n_ 

2^ 



POO 

Jo 



M-i(C)- 



(21) 



bears a remarkable similarity to ( |20| i, but now with a different 
matrix M entering into the matrix manipulations; whereas 
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Fig. 6. Repulsive Casimir force between metallic objects in vacuum. Plotted 
is the z-directed force on an elongated nanoparticle above a circular aperture in 
a metallic plate (inset), as a function of the separation distance d between the 
center of the nanoparticle and the center of the plate. The dashed red curves 
are for the case of perfectly conducting materials (for two different plate 
thicknesses), while the solid blue curve is for the case of finite-conductivity 
gold. The shaded region of the force curve indicates the repulsive regime, in 
which the nanoparticle is repelled from the plate. (The dashed vertical line 
denotes the separating plane, i.e. the value of d beyond which the nanoparticle 
is entirely above the plate.) For comparison, the dashed grey curve indicates 
the force on a perfectly-conducting spherical nanoparticle; in this case the 
force is attractive at all separations. (Figure reproduced from Ref. |11|.) 



M in (20) describes the interactions between incoming and 
outgoing waves in a muUipole expansion of the electro- 
magnetic field, M in ^2A\ describes the interactions among 



surface currents flowing on the surfaces of the interacting 



objects in a Casimir geometry. [M(^) in (21 ) is just the usual 
impedance matrix that enters into the PMCHW formulation 
of the boundary-element method [ [69| , but now evaluated at 
imaginary frequencies.] 

As one example of the type of calculation that is facilitated 
by FSC Casimir techniques, Ref. [ [TT| investigated the Casimir 
force on an elongated nanoparticle above a circular aperture 
in a metallic plate and identified a region of the force curve 
in which the force on the particle is repulsive (Fig. [6]). This 
geometry is notable as the only known configuration exhibiting 
repulsive Casimir forces between non-interleaved metallic ob- 
jects in vacuum. (On the other hand, repulsive forces between 
dielectric objects immersed in a dielectric liquid have long 
been known to exist and were observed experimentally in 
2009 1 70 1 ; in addition, numerical Casimir tools have been used 
to demonstrate theoretically the possibility of achieving stable 
suspension of objects in fluids | [7T| , and further work in this 
area may have applications in microfluidics.) 



V. Summary and Outlook 

Despite spending most of its history confined to the realm 
of pure physics, the theory and experimental characterization 
of fluctuation-induced electromagnetic phenomena is at last 
poised to take on a new role as a growth area in electrical 
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engineering. The growing ease and ubiquity of nanotechnology 
are making near-field radiative transport and Casimir forces in- 
creasingly relevant to the technologies of today and tomorrow, 
with a corresponding imminent need for engineers to account 
for these phenomena in their designs. In this connection it 
is convenient that a host of powerful computational methods, 
inspired by techniques of classical computational electromag- 
netism but extending these methods in several ways, have been 
developed over the past several years to model various fluctua- 
tion phenomena. We hope to have convinced the reader that the 
sudden conjunction of new theoretical techniques, increasing 
experimental relevance, and the paucity of known results have 
created burgeoning opportunities for computational science — 
indeed, in fields where two spheres represent a novel geometry, 
the untapped frontiers of design are vast and inviting. 

What lies in store for the future of this field? The work 
reviewed in this article has answered many questions, only to 
pose many more to be addressed in the coming years. Here we 
give a brief flavor of some challenges that lie on the horizon. 

General-basis trace formulas for heat transfer. Unlike 
Casimir forces, the theory of near-field radiation does not 
yet benefit from a compact trace formula that applies to an 
arbitrary localized basis. Existing approaches either require 
the intermediary of a spectral incoming/outgoing wave basis 
(such as cylindrical or spherical waves) that may be ill-suited 
for irregular geometries, or large-scale computations involving 
costly integral evaluations. Is a synthesis (in the spirit of the 



FSC approach of Section [IVT| ) possible or practical, and what 
form does it take? 

Fast solvers. To date, practical applications of integral- 
equation Casimir techniques have evaluated the matrix oper- 
ations in equation pT] ) (matrix inverse, matrix multiplication, 
and matrix trace) using methods of dense-direct linear alge- 
bra. These methods are appropriate for matrices of moderate 
dimension {D ~ 10^ or less), but for larger problems the 
0{D'^) memory scaling and 0{D^) CPU-time scaling of 
dense-direct linear algebra renders calculations intractable. 
A similar bottleneck was encountered many years ago in 
the computational electromagnetism community, where it was 
remedied by the advent of fast solvers — techniques such as 
the fast multipole |72| and precorrected-FFT |73| methods 
that employ matrix-sparsification techniques to reduce the 
asymptotic complexity scaling of matrix operations to more 
manageable levels; 0{D^/'^ log D) [74] or 0{D\ogD) \75j, 
176 1 are typical. Although such methods could, in principle, be 
applied to stress-tensor Casimir computations | [46| , ||64j, can 
they be made practical? Can they be applied to the FSC trace- 
formula approach, and with what performance implications? 

New experimental geometries. Until recently, theoretical 
techniques in fluctuation-induced phenomena lagged behind 
the forefront of experimental progress (indeed, as we have 
seen, it is only in the past few years that complete theoretical 
solutions for the simple sphere-plate geometry commonly 
seen in experiments have become available). This situation 
has recently begun to change; with a host of new computa- 
tional methods for near-field radiative transfer and Casimir 
phenomena becoming available in the past five years, we are 
entering an era in which theoretical predictions can be used 



to guide the design of future experiments — and, ultimately, 
future technologies. Such a reversal is not without precedent 
in the history of electrical engineering. Indeed, whereas the 
first computational algorithms for modeling antennas and 
transistor circuits were validated by checking that they cor- 
rectly reproduced the behavior of existing laboratory systems, 
today it would be unthinkable to fabricate a patch antenna 
or an integrated operational amplifier without first carefully 
vetting the design using CAD tools. Will the development of 
sophisticated modeling tools for near-field radiative transfer 
and Casimir phenomena transform those fields as thoroughly 
as SPICE and its descendants transformed circuit engineering? 
In the former case, can we use modeling tools to design 
efficient tip-surface geometries for thermal lithography, or to 
invent new solar-cell configurations that exploit the interplay 
of material and geometric properties to optimize power absorp- 
tion and retention at solar wavelengths? In the latter case, can 
we use computational tools to understand parasitic Casimir 
interactions among moving parts in MEMS devices — or to 
invent new MEMS devices that exploit Casimir forces and 
torques to useful ends? 

All of these are questions for the future of fluctuation- 
induced phenomena. We hope in this review to have piqued 
the curiosity of electrical engineers in this rapidly developing 
field — and to have encouraged readers to stay tuned for future 
developments. 

In closing, we note that all of the computational results 
presented in this review were obtained using freely-available 
open-source software packages for computational electromag- 
netism: MEEP, a finite-difference solver, and SCUFF-em, a 
boundary-element solver. (Both packages are available for 
download at http://ab-initio.mit.edu/wiki.) In 
addition to their general applicability to scattering calculations 
and other problems in computational electromagnetism, these 
codes offer specialized modules implementing algorithms dis- 
cussed in this article for numerical modeling of fluctuation- 
induced phenomena. 
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